The Dynamics of Sinking Satellites Around Disk Galaxies: A Poor 
Man's Alternative to High-Resolution Numerical Simulations 
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ABSTRACT 

We have developed a simple yet surprisingly accurate analytic scheme for tracking the dy- 
namical evolution of substructure within larger dark halos. The scheme incorporates the effects 
of dynamical friction, tidal mass loss and tidal heating via physically motivated approximations. 
Using our scheme, we can predict the orbital evolution and mass-loss history of individual 
subhalos in detail. We are also able to determine the impact and importance of the different 
physical processes on the dynamical evolution of the subhalos. To test and calibrate this model, 
we compare it with a set of recent high-resolution numerical simulations of mergers between 
galaxies and small companions. We find that we can reproduce the orbits and mass- loss rates 
seen in all of these simulations with considerable accuracy, using a single set of values for the 
three free parameters in our model. Computationally, our scheme is more than 1000 times 
faster than the simplest of the high-resolution numerical simulations. This means that we can 
carry out detailed and statistically meaningful investigations into the characteristics of the 
subhalo population in different cosmologies, the stripping and disruption of the subhalos, and 
the interactions of the subhalos with other dynamical structures such as a thin disk. This last 
point is of particular interest given the ubiquity of minor mergers in hierarchical models. In this 
regard, our method's simplicity and speed makes it particularly attractive for incorporation into 
semi-analytic models of galaxy formation. 



Subject headings: cosmology: dark matter - 
— methods: numerical 

Introduction 



galaxies: kinematics and dynamics — galaxies: interactions 



Over the past two decades, observational and 
theoretical progress have given rise to an increas- 
ingly detailed picture of how structure forms and 
evolves on galactic scales. The currently favored 
theoretical models are based on the concept of 
hierarchical clustering, in a universe dominated 
by cold dark matter (CDM). Galaxies are embed- 
ded within extended halos of dark matter, which 
form through gravitationally induced mergers of 
smaller-scale structure. The evolution of individ- 
ual dark matter halos and the formation of galac- 
tic structure, if any, inside the halos is strongly 
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dependent on the non-linear dynamics of gravita- 
tional collapse, the stochastic process of merging, 
and the subsequent evolution of the merged sub- 
structure. These processes are extremely challeng- 
ing to model theoretically. Consequently the exact 
role of mergers, especially minor mergers, in the 
formation and evolution of galaxies is still poorly 
understood, even though mergers are thought to 
be an ubiquitous feature of structure formation. 

The dynamics of structure formation have been 
studied extensively using numerical simulations. 
Simulations have been used to investigate the 
properties of populations of cluster, group, and 
galaxy-sized halos in their larger cosmological en- 
vironment (e.g. Jenkins et al. 1998; Jing & Suto 
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1998; Governato et al. 1998, 1999; Kauffmann ct 
al. 1999; Pearce et al. 1999; Sigad et al. 2000), the 
detailed properties of individual halos on scales 
where their substructure is resolved (e.g. Ghigna, 
ct al. 1998, 1999; Moore et al. 1999; Klypin et al. 
1999a, 1999b; Okamoto & Habe 1999; Lewis et 
al. 2000; Yoshikawa, Jing & Suto 2000; Fukushige 
& Makino 2000; Jing & Suto 2000), as well as 
the outcome of mergers between halos or sub- 
components within halos, both minor (e.g. Quinn 
& Goodman 1986; Quinn, Hernquist & Fullagar 
1993; Walker, Mihos & Hernquist 1996; Huang & 
Carlberg 1997; Velazquez & White 1999) and ma- 
jor (e.g. Barnes 1998; Naab, Burkert & Hernquist 

1999) . Existing simulations do not yet have the 
dynamic range to explore these different scales si- 
multaneously, however, leading to a fragmented 
treatment of the subject. Moreover, numerical 
simulations suffer from several other disadvan- 
tages: (1) they are very expensive computation- 
ally; (2) the detailed statistical properties of sub- 
structures within a halo may depend sensitively on 
the scheme used to identify them; (3) the evolution 
of these objects may be influenced in complex ways 
by the number of particles used to resolve them, as 
well as other numerical effects such as finite force 
resolution (e.g. Ghigna et al. 1999; Knebe et al. 

2000) . 

The above limitations are particularly frustrat- 
ing, given the recent suggestions that hierarchical 
models may be failing to match observations on 
galactic scales. One apparent discrepancy is in the 
number of satellites expected to have survived in 
a Milky Way-sized halo versus the observed num- 
ber of satellites around the two large galaxies of 
the Local Group (Klypin et al. 1999b; Moore et 
al. 1999; Bullock, Kravtsov & Weinberg 2000). 
This excess structure may be implicated in several 
other problems, including the small disk sizes pro- 
duced in hydrodynamic simulations (e.g. Navarro 
& Steinmetz 2000), and the question of disk sur- 
vival against heating in minor mergers (Toth & 
Ostriker 1992; Kauffmann & White 1993; Lacey & 
Cole 1993; Navarro, Frcnk & White 1994; Moore 
ct al. 1999). To resolve these issues, it is necessary 
to consider separately the effects of background 
cosmology, the power spectrum of density fluctu- 
ations, and most importantly the nature and dy- 
namics of substructure within galactic halos. Ac- 
complishing this goal numerically would entail ex- 



ploring a large parameter space, via ultra-high res- 
olution simulations. This proposition is, however, 
prohibitively expensive at present. 

An alternate approach to studying structure 
formation is to use semi-analytic (SA) methods, 
which combine analytic theory and numerical re- 
sults. Semi-analytic models of galaxy formation 
(e.g. Kauffmann, White & Guiderdoni 1993; Cole 
ct al. 1994; Somerville & Kolatt 1999) generate 
random realizations of merging (or "merger his- 
tories" or "merger trees") between halos based 
on Press-Schechtcr statistics (Press & Schcchter 
1974). The formation and evolution of galactic 
structure within the halos is then governed by a 
set of prescriptions which attempt to describe the 
effects of merging, hydrodynamics, shocks, dissi- 
pation, star formation and feedback. 

Semi-analytic models of galaxy formation have 
thus been used to study the galaxy luminosity 
function, the Tully-Fisher relation, the morphology- 
density relation and other global properties of 
galaxy populations (see Somerville & Primack 
1999 for a recent review, as well as Bullock et 
al. 2000 for more recent work on galactic satel- 
lites). SA methods have the advantage of being 
extremely fast compared to fully numerical sim- 
ulations, and although the results depend on the 
prescriptions adopted (Benson et al. 1999), these 
are, at least in principle, transparent, so that it is 
easy to test the consequences of modifying them. 
SA methods are therefore useful for exploring the 
relative importance of the various ingredients of 
the galaxy formation model — the background 
cosmology and the shape of power spectrum, the 
density profiles of dark matter halos, the non- 
linear dynamics of merging, the gas and radiative 
physics, star formation and feedback algorithms 
- in determining the appearance of present-day 
galaxies. 

Until recently, however, SA models of galaxy 
formation have focused on star formation and dis- 
sipative processes, and have included only a sim- 
plified description of merging, or even ignored the 
impact of merging altogether. In Diafcrio et al. 
(2000), for example, nearly equal-mass mergers 
result in the destruction of the disk, but the im- 
pact of uneven-mass mergers on the galaxy em- 
bedded in the more massive halo is ignored. This 
makes it difficult to ascertain whether, for exam- 
ple, morphology-related results are real or merely 
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artifacts of the oversimplified prescription. It 
also makes SA results hard to relate to studies 
of the dynamical evolution of individual galaxies, 
whether analytic (e.g. Dalcanton, Spergel & Som- 
mers 1997; Mo, Mao & White 1998) or numerical 
(e.g. Velazquez & White 1999). 

Studying various issues arising from hierarchi- 
cal clustering scenarios, including the distribution 
of satellites around galactic systems, the impact 
of these satellites on a thin disk, and the for- 
mation of the stellar halo from tidal debris, re- 
quires a method for determining the evolution of 
the satellites which synthesizes and generalizes the 
results of existing numerical studies, without re- 
sorting to expensive ultra-high resolution simula- 
tions. The method should take into account the 
internal structural properties of a satellite, its or- 
bital parameters and the details of its interaction 
with the main galaxy. 

We have developed a simple analytic scheme 
that addresses this need and complements exist- 
ing semi-analytic and numerical models of galaxy 
formation. We consider subhalos within a galaxy 
halo individually, following their orbits and ac- 
counting for dynamical friction, mass loss and 
heating using analytic expressions. This approach 
will allow us to carry out detailed studies of the 
properties of subhalos for a variety of cosmologies 
and power spectra. Since the physical processes 
underlying our scheme are modeled explicitly, we 
can determine their relative effects on subhalo evo- 
lution directly. Finally, as our model allows us 
to generate large numbers of realizations at little 
cost, we can do all of the above in a statistically 
meaningful manner. We also note that although 
our method is described in the context of galaxy- 
sized halos, it is completely general and can also 
be used to study the evolution of substructure in 
clusters, for instance. 

In this paper, we outline our dynamical model 
for the evolution of substructure. In section 2, we 
give a brief synopsis of the previous investigations 
of the dynamics of satellites merging with a larger 
galactic system, and then describe the theory un- 
derlying our scheme for following orbital decay, 
mass loss, and the tidal disruption of subhalos. 
In section 3, we test and calibrate our model by 
comparison with recent high-resolution numerical 
simulations of sinking satellites by Velazquez & 
White (1999) (VW hereafter). In section 4, we 



explore how tidal heating, the form of the galac- 
tic potential and the subhalo mass profile affect 
mass loss and orbital decay. We summarize our 
results in section 5. In subsequent papers we will 
combine our model of dynamical evolution with 
semi-analytic merger trees to study the halo sub- 
structure and disk heating produced in a galaxy 
halo by multiple mergers with cosmologically re- 
alistic satellites on representative orbits. 

2. Dynamics of Merging Substructure 
2.1. Background 

The first detailed study of the dynamics of 
satellites evolving within a halo containing a disk 
galaxy was carried out by Quinn & Goodman 
(1986). The study was subsequently improved 
upon by Quinn et al. (1993). These authors inves- 
tigated the problem using numerical simulations 
and found that while dynamical friction could not 
account for the anisotropy in the orbits of satel- 
lites around spiral galaxies observed by Holmberg 
(1969), the decay times were indeed short for large 
satellites, as expected from analytic arguments, 
although there was some variation depending on 
the details of the orbits. Prograde orbits close to 
the plane of the disk tended to decay the fastest. 
The studies considered only satellites on initially 
circular orbits and found that the sinking of the 
satellites took place in two steps: (1) A relatively 
slow decay of the orbital radius largely dominated 
by loss of altitude with respect to the disk; (2) 
a rapid decay in the radius once the satellite is in 
the disk plane. The results also suggested that the 
massive satellites would heat the disk appreciably, 
although the effect was complex and depended on 
the orbit and the internal structure of the satel- 
lite. Finally, the authors noted that noise in the 
simulations is a significant factor and that even 
in simulations with ~ 500, 000 particles, the noise 
would make it difficult to study the dynamics of 
the system over the Hubble time. 

In the period between Quinn & Goodman 
(1986) and Quinn et al. (1993), Toth & Ostriker 
(1992) (TO hereafter) studied the effects of minor 
mergers on the structure of the Milky Way's disk 
using a semi-analytic model for the evolution of 
the satellites, and the disk heating they produced. 
The orbital energy of a satellite lost through dy- 
namical friction was assumed to go into heating 
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the disk locally. Comparing their predictions to 
the observed scale height of the disk and to the lo- 
cal value of Toomre's stability parameter Q, they 
concluded that the Milky Way could not have ac- 
creted more than 4% of its mass interior to the 
solar circle in the past 5 Gyr. 

Several more recent studies have sought to eval- 
uate the results of TO. Walker, Mihos & Hernquist 

(1996) (WMH hereafter) and Huang & Carlberg 

(1997) (HC hereafter) used numerical simulations, 
improving upon previous work by including a re- 
sponsive halo. WMH opted to use a very large 
number of particles to reduce numerical noise and 
in turn, had to start their satellites fairly close 
in, at a distance of 21kpc from the center of the 
galaxy. They found that the evolution of the 
satellite was similar to that described by Quinn 
& Goodman (1986) and Quinn ct al. (1993). As 
for the disk, they found that accretion of a satel- 
lite with 10% of the disk mass would result in a 
significant thickening of the disk at the solar ra- 
dius. HC, on the other hand, considered satel- 
lites, starting out at larger distances, with masses 
of 10-30% of the disk mass. The internal struc- 
ture of the satellites consisted of a small tightly 
bound core embedded in an extended low-density 
envelope containing most of the mass. As a result, 
the satellites tended to be disrupted by tidal forces 
before they could heat the disk significantly. For 
completeness, we note that Weinberg (1995, 1998) 
and Sellwood, Nelson & Tremaine (1999) have also 
studied the effects of satellite-disk interactions, the 
focus of these studies being on the response of the 
disk. 

The most recent numerical study of satellite- 
disk interactions is that of VW, who studied 
mergers involving several different satellites, in- 
termediate in density between those of TO and 
those of HC, on various different orbits. They 
confirmed that Chandrasekhar's formula (Chan- 
drasekhar 1943) gives a useful approximation to 
the drag force exerted by the halo on the satel- 
lite provided the Coulomb logarithm is adjusted 
separately for each orbit. VW also found that 
the response of the disk depends partly on the 
orientation of the satellite orbit, prograde encoun- 
ters tending to heat the disk preferentially, while 
retrograde encounters tend to tilt it. They con- 
cluded that TO overestimated the magnitude of 
disk heating by a factor of 2-3 overall. 



Although much progress has been made in 
studying the effect of minor mergers on galactic 
structure, it is hard to determine the cosmological 
implications of merger simulations given the lack 
of clear, consistent, and robust results across the 
different studies. This is, first and foremost, due 
to the fact that each study has considered satel- 
lites with different internal properties. Structural 
characteristics such as the satellite's density profile 
affect the rate of mass loss and therefore the evo- 
lution of the satellite's orbit, not to mention the 
response of the disk to the satellite. In addition, 
the satellites in the different studies were started 
at different radii, and early studies did not include 
the dynamical friction produced by the halo. On a 
more practical level, the simulations were subject 
to numerical effects such as finite force resolution, 
shot noise, relaxation, and artificial heating from 
interactions between particles of different masses, 
which differed from study to study. Finally, all 
studies prior to that of VW focused on satellites 
on circular or nearly circular orbits, making it 
difficult to generalize their results to satellites ac- 
creted by galaxies in self-consistent cosmological 
settings. To analyze this prior body of numeri- 
cal work, and overcome the limitations mentioned 
previously, requires an alternative method which 
follows the relevant physical processes explicitly, 
and can generate many realizations of merging at 
little computational expense. 

2.2. The Semi-Analytic Model 

To model the evolution of a subhalo as it merges 
with a larger halo, we treat it as a spherically sym- 
metric satellite, with structural properties that 
change over time. The structure of the satellite is 
fully specified at any time by its mass, the form of 
its density profile, its initial core or scale radius, a 
limiting outer radius, and the amount of heating it 
has experienced. The density profile is initially set 
to a standard form such as a King model, or any 
one of several common analytic density profiles. If 
the satellite experiences tidal heating, however, its 
density profile may change. We do not track these 
changes explicitly (although they may be calcu- 
lated from our other results), as we are only con- 
cerned with changes in the mean density of the 
satellite within its limiting radius in our descrip- 
tion of heating and mass loss. 

To determine the satellite's orbit, we ignore 
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its spatial extent to first order, and calculate the 
trajectory of a point particle with the same to- 
tal mass, moving in the gravitational potential of 
the halo-galaxy system. This approximation is en- 
tirely sufficient as long as the scale of the satellite's 
orbit is larger than the satellite itself; if the satel- 
lite falls into the central region of the potential, 
we consider it to be disrupted in any case. The 
background potential is taken to be static in the 
present study, although in general it can be al- 
lowed to vary in a self-similar way in our code. 

To account for dynamical friction, the response 
of the halo and disk to the satellite, we impose a 
drag force on the satellite which we calculate using 
Chandrasekhar's formula. As mentioned above, 
we also adjust the satellite's total mass and modify 
its internal structure in response to tidal stripping 
and tidal heating, respectively, during the course 
of its orbit. The satellite is considered disrupted 
when it has been stripped down to its core radius, 
or when it falls into the central region of the po- 
tential. Each of these processes is described sepa- 
rately below. 

2.2.1. Dynamical Friction 

Chandrasekhar (1943) showed that a massive 
particle moving through a distribution of back- 
ground particles will generate a wake. The col- 
lective gravitational force from the wake will act 
back on the massive particle, causing a drag force 
know as dynamical friction. 

Dividing the background potential into its two 
kinematically distinct components, we use Chan- 
drasekhar's formula to estimate the dynamical 
friction exerted by the halo/bulge system and by 
the disk on an orbiting satellite: 

Fdf = Fdf,halo + Fdf,disk 

= -4nG 2 Ml t 

xJ2 Pi(<^,i) In A,-^* (1) 

i=h,d I V rel,»| 

where V re i,h = V sa t, V re i,d = V sat - V ro t, 

Pi{< V IcU ) = Pi (r) [erfpfO - X % erf (JQ)] , 

and Xi = |V reM |/(\/2ai). 

Here M sat is the mass of the satellite, r is its 



position, V sat is its velocity, V rot is the local cir- 
cular velocity of the disk, p^ is the local density 
of the spherical (halo/bulge) component, pd is the 
local density of the disk, In Ah and In Ad are the 
Coulomb logarithms for the halo/bulge and the 
disk, and and ad are the one-dimensional ve- 
locity dispersions of the halo/bulge particles and 
the disk particles respectively (Binney & Tremaine 
1987). 

The derivation of this formula assumes a mas- 
sive point particle, moving through an infinite, ho- 
mogeneous background of much lighter particles 
with an isotropic Maxwellian velocity distribution 
of zero mean. Numerous detailed studies of satel- 
lite dynamics (Weinberg 1986; Cora, Muzzio & 
Vergne 1997; Bontekoe & van Albada 1987; van 
den Bosch ct al. 1998; VW; Colpi, Mayer & Gov- 
ernato 1999), have shown the formula to be more 
widely applicable, however, in the sense that it 
gives a useful approximation to the drag force on 
an extended satellite in a finite halo-galaxy sys- 
tem provided that the Coulomb logarithms are 
adjusted appropriately. In general, Binney & 
Tremaine (1987) suggests that Chandrasekhar's 
formula will be fairly accurate provided that the 
mass of the satellite does not exceed 20% of the 
mass of the larger system, and that the orbit of 
the satellite lies neither outside the larger system 
nor completely within its core. 

The argument of the Coulomb logarithm can 
be expressed as A = & max /& m i n , where 6 max and 
frmin are measures of the maximum and the min- 
imum impact parameters of the background par- 
ticles contributing to the wake. For a finite back- 
ground system, 6 max is conventionally taken to be 
the characteristic scale of the system. Possible 
choices for a spherically symmetric system include 
the half-mass radius of the system (e.g. Quinn & 
Goodman 1986), the distance over which the back- 
ground density changes by a factor of two (Bin- 
ney & Tremaine 1987), and the tidal radius of the 
halo or the distance between the satellite's posi- 
tion and the center of background system (Colpi 
et al. 1999). 

The value of 6 m i n is equally ambiguous. For 
a point mass satellite, 6 m ; n = G(M sat + m)/V 2 , 
where m is the mass of the background particles 
and V is a velocity "typical" of the encounter, such 
as the r.m.s. velocity of the background particles 
(Chandrasekhar 1943), or their mean velocity rel- 
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ative to the satellite (Binney & Tremaine (1987). 
For extended satellites, White (1976) derived an 
expression for 6 m i n which is approximately equal 
to 0.2 r t (or very roughly the half-mass radius) 
for a wide range of King profiles, while Quinn & 
Goodman (1986) take b m i n to be the larger of the 
half-mass radius of the satellite and the point- mass 
value G(M sat +m)/V 2 , with V taken as the mean 
velocity of the satellite with respect to the back- 
ground particles. 

The choice of an appropriate Coulomb loga- 
rithm to describe friction from the disk, and more 
generally the applicability of Chandrasekhar's for- 
mula to an inhomogeneous distribution of back- 
ground particles, is even less clear. Maoz (1993) 
and Dommguez-Tenreiro & Gomez-Flechoso (1998) 
derived formulae for the magnitude of the energy 
loss produced by an arbitrary distribution of uni- 
form velocity dispersion, but could not specify the 
direction of the corresponding frictional force. It is 
also possible to calculate dynamical friction for a 
uniform background of particles with an ellipsoidal 
velocity distribution (Binney & Tremaine 1987), 
in which case the frictional force is strongest in 
the direction of the smallest principal axis of the 
distribution. Since disk friction is of secondary 
importance compared with halo friction, however, 
we shall limit ourselves to using Chandrasekhar's 
formula to calculate its approximate direction and 
magnitude. 

As noted above, in calculating dynamical fric- 
tion we have and taken the density of the satel- 
lite's wake to be constant, and equal to the back- 
ground density at the centre of the satellite. Since 
the wake has a finite extent, this approximation 
may result in errors in the drag force if the back- 
ground density changes over small scales. This 
is likely to occur, for example, when the satel- 
lite is in the plane disk, because of the latter's 
small vertical scale height. To correct for this, 
the disk density pd used in equation (1) should 
be smoothed in the vertical direction. In princi- 
ple, the smoothing length ought to be related to 
the characteristic scale of the wake or the satel- 
lite; however, this would mean using a different 
smoothing scale for each individual satellite. At 
present this level of complication does not seem 
warranted, nor would it fully account for finite- 
size effects (see Dominguez-Tenreiro & Gomez- 
Flechoso 1998 for a more detailed discussion of the 



drag force on extended objects). We have there- 
fore chosen to smooth the disk density in the ver- 
tical direction by a fixed length corresponding to 
two times the disk scale height, noting that this 
smoothing length is on the order of the half-mass 
radius of satellites with masses in the range where 
dynamical friction has a substantial effect. 

Given the uncertainties associated with calcu- 
lating the value of the Coulomb logarithm, on the 
one hand, and the fact that, on the other hand, 
Chandrasekhar's formula with an appropriately 
adjusted Coulomb logarithm gives an excellent ap- 
proximation to the drag force seen in numerical 
studies, we will treat In Ah and In Ad as free pa- 
rameters. We adjust their values in order to lo- 
cate a point in the In Ah- In Ad plane for which 
our semi-analytic results provide the best overall 
match to a series of fifteen satellite simulations 
carried out by VW, the most detailed of such simu- 
lations carried out to date. One could in principle 
tune the logarithms to fit each simulation sepa- 
rately. Instead, we prefer to identify a single set 
of values that works well for all the orbits, since 
we would like to determine values of In Ah and 
In Ad that can be used in more general studies of 
satellite-disk interactions. 

2.2.2. Mass Loss 

The magnitude of the drag force on a satellite 
due to dynamical friction depends on its mass. A 
satellite with a finite extent will lose mass as it or- 
bits within the halo-galaxy system, and as it loses 
mass, the drag force it experiences will decrease. 
As a result, mass loss can significantly alter the 
dynamics of the satellite. To account for this, we 
need to estimate the amount of the mass that re- 
mains bound to the satellite throughout the course 
of its orbit. 

Material becomes unbound from the satellite 
through the action of tidal forces. Slowly vary- 
ing and rapidly varying tidal forces will affect the 
satellite differently. In a slowly varying system, 
material outside some limiting, "tidal" , radius will 
be stripped from the satellite, while in a rapidly 
varying system, material throughout the satellite 
will be tidally heated. These two regimes have 
been studied previously by making the approxi- 
mation that the system is static in the first case 
(that is a satellite on a circular orbit), or that the 
satellite undergoes a very short perturbation but is 
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otherwise isolated in the second case (the impulse 
approximation). In this section, we will consider 
tidal stripping on general orbits. Tidal heating 
will be treated separately in the section that fol- 
lows. 

For a satellite on a circular orbit of radius r 
within a spherically symmetric mass distribution, 
the combined potential of the entire system is 
static in the rotating frame. In this case, we can 
identify the tidal radius with the distance to the 
saddle point in the potential interior to the satel- 
lite's orbit, since this is the point where the radial 
forces on a test particle cancel out (von Hoerner 
1957; King 1962; Binney & Tremainc 1987). The 
distance R t from the satellite center to this point 
is: 



GM S , 



u 2 - d 2 $/dr 5 



1/3 



(2) 



(King 1962), where M sat is the mass of the satel- 
lite, co is its angular velocity and $ is the potential 
of the main system. 

This estimate of the tidal radius is only for- 
mally valid when M sat is much smaller than the 
mass of the main system, i? t is much smaller than 
the orbital radius, and the satellite is corotating at 
its orbital frequency. Even under these restricted 
assumptions, the mass inside R t is only approxi- 
mately equal to the bound mass, because there ex- 
ist orbits that extend beyond Rt but remain bound 
to the satellite (Binney & Tremaine 1987, and ref- 
erences therein). Furthermore, even in this simple 
case, the tidal boundary is not spherical and thus 
the use of expression (2) is approximate. 2 

General satellite orbits are neither circular, nor 
is the external potential in which they are moving 
necessarily spherical. We can still use equation 
(2) to define an instantaneous tidal limit for the 
system, where u> is now the instantaneous angular 
velocity of the satellite. For non-circular orbits, 
or orbits out of the plane of the disk, R t changes 
with time, however, and mass outside R t will be- 
come unbound as a result of successive accelera- 
tions over the course of the orbit, rather than be- 
ing stripped immediately. While equation (2) rep- 
resents a steady-state solution to mass loss, the 



2 Innancn, Harris & Webbink (1983), for instance, calculate 
a slightly different value for the limiting radius based on 
the length of the short axis of a tidally distorted satellite. 



characteristic timescale for transient changes in 
mass on general orbits should be the orbital pe- 
riod. To model this type of mass loss, we therefore 
assume that the satellite mass beyond Rt is lost 
over the course of one orbital period, and scale the 
mass loss in each timestep accordingly. 

In calculating d 2 $/dr 2 , we will average over the 
asphericity of the potential due to the disk com- 
ponent. We set: 



d 2 $ 
dr 2 



d 2 $ 



sph 



dr 2 



_d 

dr 



-GM(< r) 



(3) 



where $ sp h is the potential produced by a spher- 
ically symmetric distribution with mass M(< r) 
interior to r. This will be very close to the radial 
gradient of the actual force on the satellite when 
it is far from the disk, or when it is in the plane 
of the disk. Only when the satellite is close to the 
disk, but on an inclined orbit, will the true gra- 
dient differ substantially from this value, and in 
practice we expect tidal shocking to dominate the 
physics of the mass loss in these cases. 

We can write the stripping condition in terms 
of densities; the tidal limit occurs at the radius 
R t within which the mean density satellite ~p sat 
exceeds the density of the galaxy interior to its 
position, /5 ga i, by a factor ij: 



with 



Psat(< Rt) = VP g A< r ) 



Psat(< Rt) 
Pgal(< r ) 

r 3 ' 
GM{< r). 



(4) 
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dr 2 
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j 2 u 2 dr 2 



(5) 



where co is the instantaneous angular velocity of 
the satellite, and co c is the angular velocity of a 
circular orbit of radius r. 

This leads to a particularly simple algorithm for 
stripping satellites. First, we divide the satellite's 
orbital path into discrete sections corresponding 
to fixed timesteps. In each timestep, we deter- 
mine the tidal radius of the satellite using equa- 
tion (4). Of the material outside this radius, we re- 
move a fraction At/t orb , where At is the length of 
the timestep and t or b — 2tt/uj is the orbital period, 
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which we assume to be the characteristic timescale 
for mass loss. Finally, we treat the satellite as dis- 
rupted and set its bound mass to zero when the 
tidal radius becomes smaller than the core radius 
of the profile, although by this point it has nor- 
mally lost so much mass that the exact disruption 
criterion is unimportant in practice. As mentioned 
above, we also treat satellites which have fallen 
into the core of the bulge as disrupted, to avoid 
instabilities in the orbital calculation. 

2.2.3. Tidal Heating 

Whereas a steady or slowly varying tidal field 
will result in the stripping of loosely bound mass, 
a rapidly changing gravitational field, caused, for 
example, by fast encounters with the galactic disk 
or bulge, will induce gravitational shocks that can 
add energy to the satellite, changing its structure 
and accelerating mass loss (Ostriker, Spitzer & 
Chevalier 1972; Spitzer 1987; Kundic & Ostriker 
1995; Gnedin & Ostriker 1997, 1999; Gnedin, 
Hernquist & Ostriker 1999). 

Tidal heating from shocks changes both the 
mean and the dispersion of particle energies within 
the satellite; to model heating fully requires a 
Fokker-Planck code to track the changing dis- 
tribution function. In keeping with the simple 
method for estimating tidal mass loss developed 
above, we will derive a first-order correction for 
tidal heating, and scale it to match the mass 
loss rates seen in the simulations. To do this we 
first identify rapid shocks by comparing the shock 
timescale to the satellite's internal orbital pe- 
riod. Specifically, we heat the satellite only when 

ishock < £orb,sat, where tshock = (* s h!d + ^sh!b) _1 

is an average of the disk and bulge shock times 
i s h,d = z/V ZiBB , t and t sh ,b = r/V sat , weighted so 
that the shorter time dominates, and £ rb.sat = 
27rrh/V c (rh) is the orbital period of the satellite 
at its half mass radius. This corresponds to the 
range of shock timescales considered by Gnedin 
& Ostriker (1999). Over the course of each rapid 
shock, we then calculate the first-order change in 
energy within the satellite, and estimate how the 
satellite's density profile will change as a result of 
this energy input. 

In order to model the effect of tidal heating on 
the satellite, we use the formalism of Gnedin et al. 
(1999). Consider an element of unit mass, with 
coordinates x with respect to the satellite center. 



In the impulse approximation, tidal acceleration 
acting over the course of a rapid encounter, of du- 
ration t, will induce a velocity change: 



AV = / A tid (t')dt' 
Jo 



(6) 



relative to the satellite's center of mass, where A t id 
is the tidal acceleration. The resulting first-order 
change in its energy is simply equal to the work 
done by the tidal forces: 

AE^t) = W tid (t) = ^AV 2 

= ± f AUt'W-f A tid (t")dt".(7) 
1 Jo Jo 

If we divide the shock into a series of n discrete 
time steps of length At, then the work done is: 



W tid (t n ) = -At 2 



^A tid (^)-^ Atidfo) 

3=0 



i=0 



■ (8) 



In going from t n to t n+ \, the energy change in a 
single timestep is: 



AWtid(t n -» t„+i) 
= ^At 2 A tid (t n )- 



2^2 Atid(ii) + Atid(tn) 



7=0 



■ (9) 



If the satellite is sufficiently small, we can ex- 
press the tidal force in terms of the gradient of 
the gravitational force due to the external poten- 
tial, evaluated at at the center of the satellite: 

Atid(t) = x(t) • [Vg]( x=0 ) = g a ,bXb{t) e a , (10) 

where g is the external gravitational field, g a .b = 
dg a /dxb evaluated at x = 0, e a is the unit vector 
in the x a -dircction and repeated indices a, b indi- 
cate summation over the three Cartesian coordi- 
nates. Thus taking the dot product in equation 
(9) and averaging over a sphere of radius r gives: 



AW ti d(t n -> Wi) 
1 



= -r 2 At 2 
6 



2g a ,b{t n ) ^ 9a,b{U) 



i=0 



9a,b(tn)g a ,b(tn) 



(11) 
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with eighteen terms from the summation over a 
and b, where we have used the fact that 

(xiXj) = —r Sij 

averaged over a sphere. (As explained below, this 
is only strictly true if the shock is rapid, so that 
Xi(to) — Xi(tj) ~ Xi(t n ) for all tj < t n .) 

There arc two important corrections to equa- 
tion (11). First, our calculation is based on the 
impulse approximation, that is the mass element 
is assumed to remain stationary over the course 
of the shock. This approximation is expected to 
break down in the central regions of the satellite 
where the dynamical timescales can be compara- 
ble to or shorter than the shock duration. In these 
regions the effects of the shock will be greatly re- 
duced. To account for this, we adjust the heat- 
ing during rapid shocks using the first-order adi- 
abatic correction discussed by Gnedin & Ostriker 
(Gnedin & Ostriker 1999, and references therein): 

AE\ = Ai(a;)A£i i i m p , (12) 

where A\(x) = (1 + x 2 )~ 7 and x, the adiabatic 
parameter, is the ratio of the shock duration i s hock 
and the orbital period of the satellite at its half- 
mass radius t rb,sat- Since most of the heating in 
our model comes from fairly rapid disk shocks, wc 
use a value of 5/2 for 7, in keeping with the results 
of Gnedin & Ostriker (1999). 

Second, heating also leads to a change in the 
internal velocity dispersion of the satellite, as dis- 
cussed by Kundic & Ostriker (1995). Both the 
overall energy gain and the increase in the dis- 
persion will cause some of the mass to become un- 
bound. In keeping with the simplicity of our semi- 
analytic approach, we only compute the first-order 
gain in energy but account for the higher-order ef- 
fects through the introduction of a heating coef- 
ficient, eh, that we will adjust to yield reasonable 
overall matches to the VW simulations: 

AE = e h A£i (13) 

where 

A-Ei = Ai(x)AEi,i mp = ii(i)AW tid . 

Kundic & Ostriker (1995) estimate that the 
second-order heating term has an effect compa- 
rable to or greater than that of the first-order 



term, so we expect eh to be greater than 2. From 
the disruption timescale arguments in Gnedin & 
Ostriker (1997), for instance, we might expect that 
eh — 7/3. The value of eh used in practice, how- 
ever, will also depend on the shocking criterion 
and the adiabatic parameters, as explained below. 

To determine how heating affects the satellite, 
we assume that the change in its mass distribu- 
tion does not involve shell crossings and that the 
potential energy of a mass element remains pro- 
portional to its total energy (as it would in virial 
equilibrium). The total energy E(r) of a mass el- 
ement at a radius r will thus be proportional to 
— 1/r, so that an injection of energy AE(r) will 
result in a change in the radius Ar cx AE(r) r 2 . 
In the absence of shell crossings, the mean density 
inside radius r will therefore change as 

. /3M(< r)\ Ar AE(r) 
*Pr= A (-7^J oc-^a— Ai. 

(14) 

As equation (14) suggests, heating will cause 
the mass distribution to expand. Some of the 
material near the tidal radius will therefore cross 
this boundary and may be stripped away. Conse- 
quently, heating will accelerate mass loss. Since 
AWtid(r-) / Y 2 is independent of r, if we keep a run- 
ning total of the eighteen terms in equation (11), 
we can calculate the approximate density change 
produced by tidal heating at some arbitrary radius 
r as a function of time, and then apply tidal strip- 
ping (eq. [4]) to the new, heated density profile to 
determine how much mass is lost. 

This method for describing heating suffers from 
two limitations. First, we have used an average 
adiabatic correction for the system in equation 
(12). The actual correction for an orbit of radius 
r will depend on the orbital period at that radius, 
so the density change produced by heating will 
also depend weakly on radius. If we use a single 
scalar quantity to track AE/r 2 for a given satel- 
lite, we will overestimate the heating experienced 
in its inner regions as a result. Second, we have 
assumed that the internal structure of the satellite 
doesn't change in the derivation of equation (14). 
On slow orbits, satellite structure may be partially 
re-virialized as the system relaxes between shocks, 
producing a tightly bound core which is resistant 
to subsequent tidal effects. In practice, we ex- 
pect these effects to be secondary, and to be partly 
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masked by uncertainties in our choice of values for 
In A and eh- 

3. Comparison with Numerical Results 
3.1. Simulation Parameters 

To set the values of the three free parameters 
(In Ah, In Ad and eh) in our semi-analytic scheme 
and to test how well this prescription does in pre- 
dicting the evolution of a subhalo moving inside 
a larger halo, we compare our model results to 15 
recent high-resolution simulations by VW of the 
evolution of a single satellite within a larger halo 
containing a disk galaxy. Reproducing the results 
of these simulations offers a good test of our sim- 
plified description of merging, since each one fol- 
lows the orbital evolution and mass loss history 
of the satellite in detail, and together they cover 
a range of different orbits and satellite densities. 
The satellites also have large masses and small or- 
bital pericenters; thus if we can match these sim- 
ulations reasonably well, we expect the agreement 
to be even better for the more common case of 
small satellites orbiting at large distances in the 
halo. 

For the purpose of comparison, we evolved or- 
bits in a static potential identical to the one 
adopted by VW, which consists of three compo- 
nents, a truncated isothermal halo with a core, a 
stellar bulge, and an exponential disk. The density 
profiles of the three components are: 

n M M H a cxp(-r 2 /rg ut ) 

- 27r3/ 2 r cut r 2 + 7 2 

, x Mb a_ 

PB[T) ~ 2tt r(a + r) 3 

Pd ( r ) = ~r cxp(- J R/i? D )sech 2 (z/z ) 
47rK D zo 

where the masses and scale lengths of the compo- 
nents are: 

M H = 7.84 x 10 n M Q ,7 = 3.5kpc,r cut = 84kpc, 

M B = 1.87 x 10 10 Mq, a = 525 pc, 

M B = 5.6 x 1O 1o M ,]?d = 3.5 kpc, z = 700 pc . 

The disk density used in the calculation of dy- 
namical friction was smoothed in the vertical di- 
rection by two disk scale heights, as explained in 



section (2.2.1), to reflect the finite size of the satel- 
lite, and we similarly smoothed the vertical com- 
ponent of the tidal field used in calculating heating 
by the disk. The sum of the halo and bulge densi- 
ties was used to calculate the other friction term, 
since these components are kinematically similar. 
For the velocity dispersions, we used: 

= (K 2 h + K 2 b ) 1/2 /V2, 

and 

ad = V Cl d/V2 = (T exp(-R/Ro) 

where V Ci h , V^b and V c ,d are the circular velocities 
of the halo, bulge and disk respectively, and we set 
a a to 143 km s -1 and R Q to 7 kpc (or 2i? d ), based 
on the velocity dispersion of the disk measured by 
VW. 

Fifteen different orbits were simulated, with ini- 
tial conditions corresponding to those of VW (see 
table 1). Our satellite models SI, S2, and S3 were 
also identical to those used by VW. These are King 
models, with core radii and initial concentrations 
appropriate to dwarf spheroidals (see table 2). 

3.2. Results 

Figure 1 shows the evolution of satellite SI on 
five orbits of different inclination with respect to 
the disk. The points are the results of the sim- 
ulation, and the solid curves are the results from 
our semi-analytic model. Figure 2 shows similar 
results for the more concentrated satellite S2. In 
each figure, the left-hand plots show the position 
of the satellite versus time, while the right-hand 
plots show the mass. The angle i indicated on the 
plots is the angle between the initial angular mo- 
mentum vectors of the satellite and of the disk, 
so that orbits G1S2 and G1S9 are coplanar with 
the disk, and prograde with respect to disk rota- 
tion, while orbits G1S6 and G1S13 are coplanar 
and retrograde. 

The semi-analytic orbits were calculated using 
Coulomb logarithms of In Ah = 2.4 for the halo 
and In Ad = 0.5 for the disk, which are in the range 
predicted by the theoretical estimates mentioned 
in section 2 (lnA h = 1.9-2.6 and lnA d = 0.6-1.3 
or less, depending on the orbit and satellite prop- 
erties). The heating coefficient used, eh = 3.0, is 
also in the expected range. We see that for this 
choice of parameter values, we obtain a very good 



10 
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Fig. 1. — Orbits and mass loss histories for satel- 
lite SI (lines), compared with numerical results 
from Velazquez & White (1999) (points). The pa- 
rameter values used were (In Ah, In Ad, £h) = ( 2.4, 
0.5, 3.0). 
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Fig. 2. — As figure 1, but for a more concentrated 
satellite (Velazquez & White model S2). 

match overall to the orbital decay and mass loss 
in all ten cases. 

Examining the orbital evolution in detail, we 
note that the semi-analytic model matches the nu- 



merical results remarkably well, especially given 
that we are using a single set of parameter val- 
ues to fit results for three different satellites and 
eight different sets of initial conditions. Our pre- 
scription for dynamical friction reproduces the de- 
cay in the amplitude and period of the orbit, and 
the semi-analytic orbit remains in phase with the 
numerical results for as long as the mass loss is 
well matched, typically five or six orbital periods. 
Varying the Coulomb logarithm by 10-20% would 
produce a better match to some orbits, such as 
the retrograde, coplanar orbit G1S6, but as men- 
tioned previously, we prefer to find a single set of 
values which fit all fifteen orbits reasonably well. 
Varying the parameters by less than 10% does not 
affect the results substantially. 

Comparing the mass loss rates, we see that the 
semi-analytic model gives an excellent estimate of 
the timescale for mass loss, and predicts the bound 
mass in the simulations to within 20%, up to the 
point where the satellite has lost most of its mass. 
Our model also reproduces the dependence of mass 
loss on the orientation of the orbit for prograde 
and retrograde orbits in the disk, predicting faster 
mass loss on prograde orbits. This appears to be 
mainly the result of the stronger dynamical fric- 
tion experienced by satellites in this case. Orbits 
out of the plane of the disk show a weak depen- 
dence on inclination in the simulations. We re- 
produce this marginally, though the amplitude of 
the effect is much smaller in our model than in 
the simulations. This may indicate that dynami- 
cal friction from the disk is more important than 
we predict for these orbits. 

Apart from considering a single satellite on a set 
of similar orbits at different inclinations, we can 
also consider different satellites on similar orbits 
or the same satellite on different orbits. Figures 
3 and 4 show the results for three different satel- 
lites: the fiducial satellite (SI), a satellite which 
is more concentrated (S2), and one which is more 
massive and more concentrated (S3), on prograde 
(Fig 3) and retrograde orbits (Fig 4), inclined by 
45°. For both prograde and retrograde orbits, the 
more massive satellite experiences more dynami- 
cal friction, falls in faster, and is disrupted. The 
more concentrated satellite retains its mass longer 
than SI, despite having fallen further into the po- 
tential. The semi-analytic model accurately repro- 
duces these trends, although for the more concen- 
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trated satellite the mass loss rates are a bit slow. 
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Fig. 3. — Orbits and mass loss histories for three 
different satellite models with the same initial or- 
bital parameters. 
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Fig. 4. — As figure 3, for a retrograde orbit. 

Dynamical friction, tidal limits, heating, and 
mass loss timescales will also depend on the circu- 
larity of a satellite's orbit. Figure 5 shows results 
for the same satellite model (SI), on inclined or- 
bits of three different eccentricities. Here again, 



we achieve an excellent match to the simulation 
results, reproducing the trend of faster mass loss 
for more eccentric orbits, and getting a good esti- 
mate of the disruption times for the three different 
orbits. 
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Fig. 5. — Orbits and mass loss histories for satel- 
lite SI, on orbits of different circularity. 

Some systematic differences are apparent in 
the comparison between the numerical and semi- 
analytic results. The mass loss in the simulations 
is a bit smoother than in the semi-analytic model, 
showing less variation in rate at pericentric pas- 
sage. We also underestimate the mass loss rates 
for the more concentrated satellite, when it is on 
orbits inclined with respect to the disk (G1S10- 
G1S12). Varying the parameter values suggests 
that this may be due to a slight underestimate 
of the dynamical friction in these cases, although 
our model for mass loss is no doubt partly respon- 
sible for the mismatch. Some of the theoretical 
estimates for dynamical friction mentioned in sec- 
tion 2 do predict a larger Coulomb logarithm for 
more concentrated satellites, so we may be seeing 
evidence of this in our results. In the absence of 
more numerical results to confirm this dependence 
on concentration, however, we will limit ourselves 
to a single set of values for the Coulomb loga- 
rithms. Finally, our most circular orbit (G1S8) 
experiences slightly less mass loss than predicted. 
In this case, the characteristic timescales for the 
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shocks are very close to the internal orbital pe- 
riod of the satellite. Using a more restrictive def- 
inition of rapid shocks in this case produces re- 
sults which match the numerical behaviour ex- 
actly. Here again though, we have insufficient nu- 
merical data to justify a general modification to 
our scheme. 

We also note that the details of the numeri- 
cal mass loss histories may depend on the precise 
definition of the bound mass used in interpreting 
the simulations. Investigating mass loss in detail 
would require a careful re-analysis of the simu- 
lations, given the importance of this and other 
purely numerical effects on the mass loss histo- 
ries. VW show, for instance, that changing the 
resolution of a simulation by a factor of four can 
have an effect comparable to the discrepancy we 
see between the semi- analytic and numerical re- 
sults (VW, figure 10). The softening lengths and 
time-stepping algorithms used could affect the nu- 
merical results to a similar degree. It is also in- 
trinsically hard to separate out the effects of tidal 
heating and tidal stripping in the results of VW, 
since the mass loss timescales produced by the two 
are similar in their simulations. The fact that we 
find good agreement with their results over a range 
of orbits and for several satellite models, however, 
gives us some confidence in our description of these 
phenomena. 

In summary, using a simple model of dynam- 
ical friction, tidal heating and tidal mass loss, 
we can reproduce analytically the results of high- 
resolution V-body simulations of mergers, includ- 
ing accurate timescales for satellite infall and dis- 
ruption, with the correct dependence on satellite 
mass and concentration, and on orbital param- 
eters. Our model has three free parameters - 
In Ah, In Ad, and eh- Of the three, the results de- 
pend most strongly on In Ah and eh, with In Ad 
making a secondary contribution. Requiring that 
our results match those of VW fixes the values of 
these parameters to within roughly 10%. The val- 
ues we obtain fall within the range predicted by 
first-order estimates of friction and heating. 

4. Discussion 

4.1. The importance of tidal heating 

Since the processes underlying dynamical evo- 
lution are specified explicitly in our model, it is 



possible to adjust them to test their relative con- 
tribution to satellite evolution. In particular, we 
can test the importance of tidal heating, often ne- 
glected in the study of satellite dynamics, on our 
results. The solid curves in figure 6 show sev- 
eral orbits calculated without tidal heating (solid 
curves), compared with the heated orbits (dotted 
curves), and the simulations (solid dots). We see 
that the overall effect of heating is to increase mass 
loss, which in turn reduces dynamical friction. In 
general, the simulation results are better matched 
by including heating, although the importance of 
heating varies with circularity and inclination. On 
inclined orbits, satellites are strongly affected by 
heating, while its effect on orbits in the plane of 
the disk is minor. This demonstrates that for the 
orbits we have considered, disk shocks dominate 
over bulge shocks as a source of heating. 
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Fig. 6. — A set of orbits calculated with and with- 
out heating (dotted and solid lines respectively), 
compared with simulations (points). The param- 
eter values used were (In Ah, In Ad, eh) = ( 2.4, 
0.5, 3.0) and ( 2.4, 0.5, 0.0) for the heated and 
unheated cases, respectively. 

If we consider a satellite to be disrupted when 
it has lost some large fraction of its mass, the dis- 
ruption times we measure for our satellites are up 
to 40% shorter due to heating. One might expect 
an even stronger effect, but we note, comparing 
the mass loss and orbital decay curves, that dy- 
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namical friction conspires to reduce the difference 
between the mass loss times for the satellites and 
orbits considered here. In the no-heating runs, 
the satellites retain more of their mass initially 
and therefore experience more dynamical friction. 
The satellite orbits decay faster as a result, and 
the satellites fall into the center of the potential 
and are disrupted. This accounts for the sharper 
cutoff to some of the mass loss curves. If the over- 
all timescales for disruption are similar with and 
without heating, this is partly an accident of the 
concentrations and masses of our satellites. A less 
massive satellite of similar density, for instance, 
would have a much longer orbital decay time but 
would lose mass through heating at about the 
same rate. More concentrated satellites may resist 
heating almost completely, while less concentrated 
satellites may quickly be disrupted by it. Finally, 
heating will produce a quite different distribution 
of stripped material from satellites, which is im- 
portant in considering halo substructure formed 
by satellite accretion and disruption. Thus, heat- 
ing cannot be neglected in studying minor mergers 
with semi-analytic models. 

4.2. The effect of the disk and the bulge 

The disk was shown above to have a strong 
effect on satellite evolution by heating satellites 
on inclined orbits. Its presence should also gener- 
ate dynamical friction, particularly for prograde, 
coplanar orbits. The effect of the bulge is less 
clear. Understanding the role of these structures is 
important in relating small-scale simulations such 
as those of VW to larger cosmological simulations, 
which do not yet include disks or bulges. To test 
the effect of these components on satellite orbits, 
mass loss, and disruption times, and to determine 
any systematic trends affecting simulation results, 
we have run our model with one or both of these 
components removed. 

Removing the bulge from the potential has lit- 
tle effect on the satellite orbits, acting only to 
decrease friction slightly. We expect the disk to 
have a much greater effect, due to its greater mass, 
which produces more dynamical friction, and to its 
vertical density gradient, about ten times larger 
than that of the bulge or halo, which should con- 
tribute roughly 100 times more heating than the 
other components, for satellites crossing the disk 
plane. Running the model without a disk confirms 



that this is indeed the case. 

In figure 7, we show satellite orbits G1S2-G1S9, 
recalculated in the same potential without a disk 
(dashed curves) , as well as the previous results for 
orbits in the presence of a disk, but with heat- 
ing turned off (dotted curves), and with both a 
disk and heating (solid curves). When the disk is 
absent, we see that the dependence of orbital evo- 
lution on inclination vanishes, as expected. Fur- 
thermore, the initial mass loss rate is reduced, and 
satellites fall further into the potential without los- 
ing as much mass. Turning off the disk or turn- 
ing off heating produces similar results for orbits 
G1S3-G1S5, indicating that the effect of the disk 
is mainly to heat satellites on inclined orbits. In 
the prograde, coplanar orbit G1S2, on the other 
hand, the disk contributes mainly to dynamical 
friction. For this orbit, turning off heating does 
not change the results substantially, whereas turn- 
ing off the disk does. 
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Fig. 7. — A set of orbits calculated without a disk 
component in the potential (dashed lines), com- 
pared with the no- heating case (dotted lines) , with 
the full model (solid lines) and with simulations 
(points). The parameter values used were (In Ah, 
In Aa, e h ) = ( 2.4, 0.5, 3.0 or 0.0). 

Overall, we conclude that the presence of disk 
has an important effect on the evolution of satel- 
lites on orbits with pericenters of 20kpc or less 
(about 6 scale radii). For the typical eccentric- 
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ities seen in cosmological simulations (Ghigna et 
al. 1998), this suggests that orbits with apocenters 
of 120 kpc or more will be affected by the disk. The 
effect of the disk on the disruption times we mea- 
sure is to reduce them by 20-30%, but as with 
heating, this difference could be much larger for 
satellites of different masses or concentrations. To 
produce realistic distributions of galactic satellites 
by semi-analytic means, (e.g. Bullock et al. 2000), 
one should therefore account for the effects of a 
disk in the model. 

4.3. Results for different satellite profiles 

Finally, the mass and concentration of a satel- 
lite can substantially affect its dynamical evolu- 
tion. One difficulty in comparing numerical stud- 
ies of disk heating through minor mergers is the 
fact that different authors have considered differ- 
ent satellite models in their simulations. In this 
section, we recalculate a few orbits using some of 
the models that appear in the literature, to test 
the effect of a satellite's internal structure on its 
orbital evolution. 

Figures 8 and 9 show several orbits from the 
set used previously, recalculated for five different 
satellite models similar to those used in recent 
merger simulations. In figure 8, the solid curves 
are for the VW satellite SI, the dotted curves are 
for a more concentrated king model, similar to the 
one used by HC, and the short-dashed curves are 
for the highly concentrated satellite considered by 
TO. In figure 9, the solid curves are for VW SI as 
before, the long-dashed curves are for the satellite 
model used by WMH, and the dot-dashed curves 
are for a satellite with an NFW profile (Navarro et 
al. 1996, 1997) of comparable concentration. All 
of the satellite models have been given the same 
mass for purposes of comparison. Their density 
profiles and structural parameters are listed in ta- 
ble 2. 

We see from the mass loss curves the strong 
effect that the satellite model has on dynamical 
evolution. The HC satellite is much more concen- 
trated than SI, but has a similar core radius, and 
by implication an extended diffuse envelope, con- 
taining most of the satellite's mass. This diffuse 
material is stripped off early on in its orbit, lead- 
ing to much slower orbital decay. The TO satel- 
lite behaves in the opposite way — its mass is so 
tightly bound that it experiences almost no tidal 
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Fig. 8. — A set of orbits calculated for various 
satellite models. The solid lines are for Velazquez 
& White's satellite SI, the dotted line is the satel- 
lite of Huang and Carlberg, and the short-dashed 
line is the more concentrated of Toth and Os- 
triker's two satellites. See text and table 2 for 
the details of the models. The parameter values 
used were (lnA h , lnA d , e h ) = ( 2.4, 0.5, 3.0), as 
above. 

stripping, falling directly into the center of the po- 
tential with little mass loss. Finally the WMH 
satellite looses a comparable amount of mass to 
SI initially, but its overall mass loss history pro- 
duces much slower orbital decay. The same is true 
for the satellite with an NFW profile. 

These results are consistent with those of the 
original studies. HC found that their satellites 
were stripped of most of their mass before they 
hit the disk, while TO's more concentrated satel- 
lite retained 90% of its mass as it fell in on a cir- 
cular orbit to a final radius of 4 kpc. WMH found 
that a fair amount of mass was stripped off in the 
outer regions of the halo-disk system, but that the 
satellite managed to carry as much as half of its 
mass into the central few kpc. Here we have com- 
pared satellites which differ only in density profile, 
on the same orbit in the same potential. Since the 
authors mentioned above consider different satel- 
lite masses, different orbits, and different forms 
of the galactic potential in their simulations, the 
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Fig. 9. — As figure 8, for Velazquez & White's 
satellite SI (solid line), the satellite of Walker, Mi- 
hos & Hernquist (long-dashed line) , and a satellite 
with an NFW profile (dot-dashed line). See table 
2 for the details of the models. The parameter val- 
ues used were (In Ah, In Ad, £h) = ( 2.4, 0.5, 3.0), 
as above. 

semi-analytic mass loss rates shown in figures 8 
and 9 will differ in detail from their numerical re- 
sults, but we certainly reproduce all of the trends 
mentioned. 

One way of understanding these different mass 
loss rates is in terms of the fraction of a satel- 
lite's mass that lies within a given mean density 
contour. This structural property is related to 
the concentration of a satellite, and to its den- 
sity profile. Figure 10 shows the mass fraction as 
a function of density, plotted for the different pro- 
files considered here. We see that all the mass 
in the HC model is at lower densities than the 
VW model SI. Almost all the mass in the TO 
profile, on the other hand, is at densities much 
higher than SI, and higher than the central density 
of the main galaxy (which is roughly 1 M Q pc -3 ). 
For the WMH profile, half the mass is at densities 
lower than VW model SI, but the core of the satel- 
lite is at higher densities. These different profiles 
lead to different mass loss rates throughout the or- 
bit of the satellite, and as a result, very different 
dynamical histories. 




6 4 2 -2 -4 -6 -8 

lo g(p„v) ( m s/p c ~ 3 ) 

Fig. 10. — A plot of the fraction of mass within a 
given radius, as a function of mean density within 
this radius, for the various satellite models. Line 
styles are as in figures 8 & 9. 

In particular, the amount of mass a satellite 
looses in the outer part of its orbit, before it hits 
the disk, can vary tremendously from one model 
to another. Figure 11 shows all the disk crossings 
recorded in three orbits (G1S1, G1S3 and G1S8), 
plotted in terms of the fraction of the satellite's 
original mass that is still bound to it at that point, 
versus the radius at which it crosses the disk. The 
different symbols indicate the satellite models of 
TO (squares), HC (triangles), and WMH (circles). 
We see that while the TO satellite crosses the 
disk many times with almost all of its mass intact, 
less dense satellites such as that of HC have been 
stripped of most of their mass after a few orbits. 
The filled symbols in figure 11 indicate the aver- 
age mass fraction for all disk crossings between 8 
and 12kpc. While TO's satellite encounters the 
disk at this radius with 96% of its mass intact on 
average, the satellites of HC and WMH have only 
20% of their mass intact at this point. Given that 
TO saw more disk heating in their study than HC 
or WMH, these results suggest that heating may 
be simply related to the mass of material accreted 
by the disk, once tidal stripping has been taken 
into account. We shall investigate this possibil- 
ity in detail in a subsequent paper. In any case, 
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when studying minor mergers, accretion and disk 
heating, it is clearly important to use cosmolog- 
ically motivated density profiles, and to consider 
how different satellite models may affect the final 
results. 
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Fig. 11.— Disk crossings in orbits G1S1, G1S3 
and G1S8. The fraction of mass remaining bound 
to the satellite is plotted versus the radius at which 
it crosses the disk. The symbols indicate the satel- 
lite models of TO (squares), HC (triangles), and 
WMH (circles). The bold symbols indicate the 
average mass remaining for all disk crossings be- 
tween 8 and 12 kpc. 

5. Conclusion 

While there has been much progress recently 
in the understanding of how structure forms in 
dark matter on cluster, group and galaxy scales, 
the dynamical evolution of halos and subhalos is 
still not very well understood. Numerical simu- 
lations typically lack the resolution and statistics 
to follow the formation and evolution of structure 
across the range of scales involved, and much of 
the underlying physics remains uncertain. Several 
observed features of galaxies, such as thin disks, 
seem difficult to explain in current hierarchical 
models. This may partly be the result of the com- 
putational limitations restricting current numer- 
ical studies, or it may imply a genuine problem 



with the underlying cosmological models. To ex- 
plore the parameter space relevant to these issues 
requires a method which is faster and less com- 
putationally intensive than numerical simulation. 
To this end, we have developed a simple model of 
the dynamical evolution of substructure on galac- 
tic scales. 

Our model follows the dynamics of individual 
subhalos numerically, but accounts for dynam- 
ical friction, tidal mass loss and tidal heating 
using analytic expressions from Chandrasekhar 
(1943), King (1962), Gnedin & Ostriker (1999) 
and Gnedin et al. (1999). We calibrate the model 
by comparison with fully numerical simulations. 
In particular, we find that we can reproduce the re- 
sults of the most recent set of high-resolution sim- 
ulations of satellite infall by Velazquez & White 
(1999), and that matching these simulations sets 
our three free parameters to within roughly 10%. 
The values we obtain are all in the range predicted 
by first-order estimates of friction and heating. 

Varying the shape of background potential, the 
amount of tidal heating and the density profile of 
the satellite in our model, we can start to extract 
from these simulations the factors contributing to 
mass loss and orbital decay. In general, tidal heat- 
ing increases the mass loss and orbital decay times 
of our satellites substantially, although these ef- 
fects are partly masked by dynamical friction for 
the satellites and orbits we consider. We find in 
particular that the presence of a thin disk strongly 
affects evolution of objects in the inner regions of 
the galaxy, while the presence of a central bulge 
has little effect. For the orbital eccentricities typ- 
ically seen in cosmological simulations, satellites 
on orbits with apocenters of 120 kpc or more will 
pass through the disk repeatedly within a Hubble 
time, so it is important to consider its effects when 
studying galactic satellites. The overall evolution 
of a satellite is sensitive to its density profile. In 
the tidal truncation approximation, for instance, 
the satellite's mass loss history is determined by 
its mass fraction as a function of mean density. 
This dependence may explain some of the discrep- 
ancies found between different simulations of disk 
heating through satellite infall. 

Given the excellent overall match to simulations 
of minor mergers that we achieve using a simple, 
physically motivated model of satellite dynamics, 
we can go on to consider the evolution of the large 
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numbers of subhalos that a galactic halo will ac- 
crete over its lifetime. In the next paper in this 
series, we describe how to construct the mass ac- 
cretion history of a halo from a merger tree, and 
use it as the input to our model of dynamical evo- 
lution. We shall then apply this method to sev- 
eral outstanding problems in galaxy formation, no- 
tably the question of disk survival in hierarchical 
models. 
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Table 1 

Summary of Simulations from Velazquez & 
White (1999) 





Name 


Satellite 8; 


f J 7\, 
C J ' p 


1 a 






model 


(kvc) 






G1S1 


SI 45° 


0.33 5.25 


59.0 




G1S2 


SI 0° 


0.55 10.5 


55.0 




G1S3 


SI 45° 


0.55 10.5 


55.0 




G1S4 


SI 90° 


0.55 10.5 


55.0 




VjrloD 


ol loo 


0.55 10.5 


55.0 




G1S6 


SI 180° 


0.55 10.5 


55.0 




G1S7 


SI 0° 


0.82 21.0 


46.5 




G1S8 


SI 45° 


0.82 21.0 


46.5 




G1S9 


S2 0° 


0.55 10.5 


55.0 




G1S10 


S2 45° 


0.55 10.5 


55.0 




G1S11 


S2 90° 


0.55 10.5 


55.0 




G1S12 


S2 135° 


0.55 10.5 


55.0 




G1S13 


S2 180° 


0.55 10.5 


55.0 




G1S14 


S3 45° 


0.55 10.5 


55.0 




G1S15 


S3 135° 


0.55 10.5 


55.0 




Note.- 


- 0i refers to the 


angle between the initial 




angular momentum vector of the satellite and the an- 




gular momentum vector of the disc. The 


circularity 




ej = J/J c , where J is the initial angular momentum 




of the orbit and J c is the angular momentum of a ci- 




cular orbit with the same energy. r„ and 


r a are the 




initial pericentric and apocentric radii of the orbit, re- 




spectively. 












Table 2 








Satellite Models 




Satellite 


Model 


Density Profile 


Mass 


r c n 








(M ) 


(pc) (kpc) 


VW SI 


King 


^/(r 2 +r c 2 ) 


5.6 x 10 9 


1000 6.3 


VW S2 


/ / 


out to r t 


5.6 x 10 9 


500 6.3 


VW S3 


/ / 


i i 


1.12 x 10 10 


850 8.5 


HC 


/ / 


i i 


5.6 x 10 9 


812 55.3 


TO 


Jaffee 


A/r 2 (r + r c ) 2 


5.6 x 10 9 


100 47.2 


WMH 


Hcrnquist 


A/r(r + r c ) 3 


5.6 x 10 9 


525 45.5 


NFW 


NFW 


A/r(r + r^) 2 


5.6 x 10 9 


500 5.0 
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